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Abstract 


The  production  of  synthetic  seismograms  using  a 
technique  originally  developed  by  Alekseev  and  Mikhail 
is  presented.  An  explicit  finite-difference  solution  t 
elastodynamic  wave  equation  for  a  vertically  inhomogen 
medium  is  determined  after  reducing  the  dimensionality 
the  wave  equation  through  the  use  of  finite  Hankel 
transforms.  The  finite  difference  scheme  used  for  both 
and  P-SV  cases  is  outlined  in  detail  as  well  as  the 
stability  criterion  for  each.  Seismic  sources  consider 
include  SH  point  torque  source,  horizontal  point  force 
vertical  point  force,  and  explosive.  The  numerical  mod 
presented  illustrate  the  nature  of  wave  propagation  fo 
particular  model  containing  two  thin  low-velocity  laye 
a  homogeneous  half-space. 
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1.  INTRODUCTION 


The  computation  of  synthetic  seismograms  using  a 
f inite-di f ference  method  is  considered  for  the  case  of  a 
vertically  inhomogeneous  medium.  Finite  Hankel  transforms 
are  used  to  reduce  the  dimensionality  of  the  wave  equation 
and  allow  the  use  of  an  explicit  finite-difference  scheme. 

1.1  Equations  of  Motion 

For  a  vertically  inhomogeneous  medium  the  equation 
describing  the  motion  is  given  by  Grant  and  West  (1965): 

1*1  pu  =  (A+y)V6  +  6VA  +  yV  u+  (VyV)u  +  V(Vyu). 

In  the  following  we  use  cylindrical  coordinates  (r,0,z) 

with  the  z-axis  pointing  downward  into  the  vertically 
inhomogeneous  half-space  whose  surface  is  at  z=0.  The 
density  is  given  by  p(z)  and  Lame's  parameters  are  given  by 
X(z)  and  y(  z).  Using  6=V*u  and  the  vector  definition  in 
cylindrical  coordinates  V2u  =  V(V*u)  -  VxVxu  this  reduces  to 

pu  =  ( A+2y ) V ( V • u)  +  ( V  »u) VA  -  yCVxVxS) 

+  (Vy  •  V)  u  +  V  (Vy  *u) 

where  "u  =  (ur,u^,uz)  and  ur  =  u^  (r ,  z  ,  t) 

u  .  =  u  (r  ,z  f  t) 

<P  <P 

U  =  u  (r ,  z  ,  t)  . 
z  z 

To  illustrate  the  decoupling  of  the  SH  and  P-SV  types 
we  may  consider  the  equation  of  motion  given  by  1.2  in 
cylindrical  coordinates: 
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V  (V-u) 


1  _3_  (ru  ) 
r  3r 
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1  9  (ru  ) 

r  8r 
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V  *  uVA 
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Figure  1.1  Coordinate  axes  orientation. 
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The  equations  of  motion  become 


Radial  Component: 


32u 


3 1' 


=  n 


3r 


(ru  ) 
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1.3 
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Vertical  Component: 


P 


(rur) 


>  +  (i  J-  (rur> 
r  3r 


+ 


3_A 
3  z 


u 

z 


Azimuthal  Component: 
1.5 


In  this  form  it  is  evident  that  SH  waves,  given  by  1.5 
separate  from  the  coupled  P-SV  waves  described  by  equations 
1  .  3  and  1.4. 


1.2  Seismic  Sources 

Body  force  terms  representing  seismic  sources  may  also 
be  written  in  terms  of  potentials  to  incorporate  in  the 
equations  of  motion.  For  SH  waves  we  may  consider  the  case 
of  a  torque  source  at  the  free  surface,  z=0.  The  body  force 
is  defined  in  terms  of  the  potential  v  as  (Aki  and  Richards 
(1980)  Vol  I ) : 
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1.6 


?  =  V  x  ( 0 , 0  ,  v )  . 


Both  Cartesian  coordinates  (x,y,z)  and  cylindrical 
coordinates  (r,0,z)  can  be  used  where  they  share  the  same 
depth  dependence: 

1.7 

v(x,y,z)  =  6  (x)  6  (y)  6  (z)  f  (t) 


v (r , 0, z) 


Itriiiz) 

2nr 


where  6  represents  the  Dirac  delta  function  and  f(t) 
represents  the  time  dependence  of  the  source.  As  a 
mathematical  model,  the  point  torque  source  is  thus 
expressed  in  terms  of  this  potential  as 
1 • 8  F(r,z,t)  =  Vx(0,0,v) 

3  v  2 

=  ’  37  * 

f  ( t)  6  ( z )  3  f  6 (r ) \ ~ 

2 it  3r  \  r  J  ® 


For  an  explosive  source,  the  body  force  term  is  defined 


as 


1.9 


F(r,z,t)  =  Vv 


=  f  (t)  V 


6 (r) 6 (z-d) 
2irr 


where  the  source  is  assumed  to  be  at  depth  z=d. 
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In  the  case  of  a  vertical  point  force,  the  force  may  be 
expressed  directly,  rather  than  in  terms  of  a  potential  as 


1.10 


F(r,z,t) 


f(t)5(r)6(z-d)  ± 
- TE -  2 


2.  SH  CASE 


2.1  Statement  of  the  Problem 

For  the  case  of  a  SH  torque  source,  only  the  u 
component  is  nonzero  due  to  the  decoupling  of  P-SV  from  SH 
motion  for  isotropic  media.  The  equation  for  u  =  u<&(r,z,t) 
becomes 


2.  1 


p (z)  32u  _  y ' (z)  3u  32u  + 

y(z)  St2  U(z)  3z  gz2 


32u  1  3u 
3r2  r  3r 


u 


subject  to  initial  conditions 


2.2 


u 


t=0 


3_u 

at 


t=0 


0. 


The  boundary  condition  imposed  by  a  torque  source  at 
the  free  surface  z=0  imposes  the  boundary  condition 
2.3 


T 


4>z 


f (t)  _d_  /  6 (r)  \ 
2  tt  dr  \  r  / 


from  which  is  obtained 
2.4 

f (t)  _d_  (  6 (r)  \ 
2tt y  dr  \  r  / 

z=0 


3u 
3  z 


The  statement  of  the  problem  for  the  case  in  which  only 
SH  waves  appear  is  given  by  equation  2.1  subject  to  initial 
conditions  2.2  and  boundary  condition  2.4.  The 
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dimensionality  of  the  equation  of  motion  can  be  reduced  by 
the  application  of  a  first  order  finite  Hankel  transform. 
The  determination  of  the  kernel  to  be  of  value  in 


transforming  the  equation  is  obtained  by  considering 
solutions  of  the  equation  Lu=0  where  L  denotes  the  linear 
differential  operator 


2.5 


Lu 


1 

p  (r) 


+  m  (r)  u. 


The  equation  of  motion  for  SH  waves  can  be  expressed  in 
terms  of  this  operator  as 


2.6 


if  ®_ 

r\3r 


V 


u  + 


_1 _  9ji  3u  _  p  (z) 

y  (z)  3z  "cTz"  y  (z) 


0. 


The  finite  transforms  corresponding  to  the  operator 


2.7 


L 


32  ,11 

T?  r  3r 


2 

v 


are  called  finite  Hankel  transforms.  If  we  consider  the 
solution  of  the  Sturm-Li ouvi lie  equation 


(K)u  =  0 


where  we  take  £=-k|2  then  the  equation  reduces  to  Bessel’s 


equat i on : 


d2u  ,  1  3u 

T- 2  +  r  3r 
9r 


2.8 


+ 


0. 


'  • 
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Imposing  the  boundary  condition  u|ra  =0  it  has  the  general 
solution  u(r)=J*>(kir)  where  k  j  are  the  roots  of  Jp  (  k  j  a ) -0 . 
( Sneddon  (1972)). 


The  transform  to  use  in  the  case  of  SH  waves,  for  v=  1  , 
is  a  finite  first  order  Hankel  transform  of  the  first  kind 


2.9 


S(ki,z,t)  = 


u (r , z , t) J1 (kir) rdr 


The  inversion  formula  is  given  by 


2.10 


00 


u  (r  ,z  ,t)  =  — 5-  l 
aZ  i=l 


S(ki,z,t) J1(kir) 


(J  2 (k±a) 


The  application  of  the  transform  is  simplified  by 
writing  equation  2.1  in  modified  form 
2.11 


u  (z) 


3 

Or 


P  (z) 


For  the  first  term,  integration  twice  by  parts  yields 
the  following  result: 


a*<z>  &(!? +  ?K(kir)rdr 


=  -p  <;  akiJQ  (k^a)  u 


r=a 


ki  s‘ki 


If  the  boundary  condition  u|rg  =  0  is  imposed,  it  is 
equivalent  to  the  introduction  of  a  totally  reflecting 
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surface  at  r=a.  If  this  reflecting  surface  is  set  far  enough 
from  the  source,  the  seismic  response  in  the  region  of 
interest  will  not  be  influenced  by  any  spurious  reflections. 
The  equation  of  motion  upon  application  of  the  transform 
becomes 


2.12 


3 

37 


y(z)  -  k2y  (z)  S  =  p(z) 


32S 

3 12 


t>0 
z>0 . 


The  boundary  and  initial  conditions  given  by  equations 
2.2  and  2.4  are  transformed  to  yield 
2.13 


3_S_ 
3 1 


t=0 


=  0 


t=0 


2.14 


3_S 
3  z 


f  (t) 


z=0 


2ny 


°  i 


3  4  i(lcir,rdr 


f  (t) 


2ttu 


o  : 


6-(— ■  k  •  J  (k  .  r)  rdr 
r  i  o  i 


f  (t)ki 
47Ty^ 


In  addition  to  the  reflections  from  the  fictitious 
surface  at  r=a,  since  the  depth  is  finite  when  the  problem 
is  solved  by  numerical  methods,  there  is  a  need  to  minimize 
reflections  from  the  lower  grid  boundary.  This  may  be 
accomplished  by  having  the  depth  great  enough  that 
reflections  do  not  influence  the  seismic  record  for  the  time 
being  considered.  Another  method  which  minimizes  such 


reflections  and  results  in  a  manageable  grid  depth  requires 
the  introduction  of  a  damping  term  7(2)  3u/9t  in  the 
equation  of  motion.  In  effect  reflections  remain  but, 
because  they  are  heavily  damped,  their  amplitude  will  be 
negligible.  The  equation  of  motion  becomes 
2.15 

h  (v(z)  H)-  kiy(z)S  =  p(z)  0+  Y(z)  tt 

where  y(z)=0  in  the  region  of  interest  and  increases 
linearly  with  depth  over  several  wavelengths  to  the  lower 
grid  boundary. 


2.2  Finite-Difference  Equation 

The  set  of  equations  to  be  solved  by  finite-difference 
methods  is  given  by 


2.15 


k^y  (2) S 


p  (z)  ^—1-  +  Y  ( 2 ) 

at 


3_S 

at 


2.13 


s 


ii 

at 


t=0 


t=0 


0 


2.14 


as 

3z 


2=0 


k±f  (t) 


where  k-,  is  the  i-th  root  of  J1(kia)  =  0.  For  a  given  root, 
this  system  can  be  solved  efficiently  by  the  use  of  an 
explicit  finite-difference  method.  It  is  necessary  to  limit 
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grid  dispersion  by  the  selection  of  a  sufficiently  fine 
spatial  grid.  The  grid  dispersion  causes  a  delaying  and 
broadening  of  the  signal  such  that  the  pulses  develop  an 
oscillatory  tail.  In  the  programs  used  it  was  necessary  to 
use  a  minimum  of  40  grid  steps  per  wavelength. 

The  determination  of  a  finite-difference  analogue  for 
equation  2.15  is  facilitated  through  the  use  of  a  relation 
given  by  Mitchell  (1977): 


2.16 


8  (  ,  .  as 


z—  z. 


*kSk-l~  (ak+ak+l>  Sk+ak+lSk+l  ,  0(Az)2 
(Az)2 


where 

2.17 


ak  =  iz 


dz 


-1 


U  (z) 


k-1 


2l^p 


k-1 
^ *  • 

k-1 


In  this  relation  the  i-th  time  step  has  been  considered  and 
the  finite-difference  formula  is  centered  about  depth  z=zk. 
The  formal  development  of  this  relationship  is  as 


a)  _  3  S 
y  ( z )  3  z 


OJ 


k-*5 


dz_  =  s  _  s 
p(z)  k  k-1 


'k-1 


follows : 


:  a 


13 


i-1  i  i  +  1 


Figure  2.1  Finite-difference  grid. 


where 


hence 


co 


k+*2 


'k+1 


dz 
y  (z) 


Sk+1  "  Sk 


a 

Tz 


(  ,  .  a  S\  _  Aco 

^(Z)  37j  -  aF 


_  1 


Z7  K+^k-^ 


”  Bk(Sk+l”Sk)  “  Ak(sk”sk-1 


B,  = 


_1_ 

Az 


'k+1 


dz 


y  (2) 


A,  = 


dz 


-1 


y  (2) 


a3z(y(z)  Tz)  Ak+1 (Sk+l”Sk}  Ak(Sk”Sk-l) 


akSk-l  ~  (ak+ak-l)Sk  +  ak+lSk+l 

(Az)  2 
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The  finite-difference  analogue  to  2.15  is  given  by 


2.18 


S1 

k+1  k+1  ?  i-l 


In  this  form  it  is  clear  that  Sk*1  can  be  solved  for 
explicitly  at  time  step  (i  +  1)  in  terms  of  known  values  of  Sk 
at  the  previous  time  steps  i  and  (i-1). 


2.19 


l<k<m 


2.3  Free  Surface  Condition 

At  the  free  surface,  k=0,  direct  application  of  2.19 
results  in  the  introduction  of  a  fictitious  grid  point 
element  Si-,.  To  overcome  this  difficulty  it  can  be  assumed 
that  the  medium  is  homogeneous  for  k=0,1,2  and  that  no 
damping  term  is  required.  With  these  restrictions  equation 
2.15  can  be  written  in  simplified  form: 


2.20 


y  (z) 
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As  a  finite-difference  equation  this  becomes 


2.21 


f  1  f 32S 

1  2  W2 

1  v  N  3t 
s 


where  v 


2  _  p (z) 
s  p  (z) 


-  2S1  + 

o  1 

(Az)2 


k2  S1 
1  o 


1 

V 


rs1+1  -  2S1  +  S1"1 
o  o  o 


t 


(At) 


To  remove  the  fictitious  value  S i  -i  boundary  condition 
2.14  may  be  used.  Its  finite-difference  equivalent  is 
2.22 


c 1  _  c1 
S1  -1 

2  Az 


f  (t)ki 
4ttu 


Substituting  for  S 1  -j  in  the  equation  of  motion  at  the  free 


surface,  the  explicit  solution  for  S 


i  +  1 


is  given  by 


2.23 


S1  +  1  =  (  2v2  (At^?  )  S1  - 
°  '  S  (Az)2  '  ° 


2  2 

2vs(At)  2  2  2 

— - 5-  +  k  v  (At)  -2 

(Azr  1  s 


-S 


i-1 


k^f  (t)  (At) 
2ttpq  (Az) 


2.4  Seismic  Pulse 


In  the  computer  programs  the  form  of  the  seismic  pulse 
f(t)  was  that  of  an  exponentially  damped  sine  function 
expressed  as 
2.24 

2TT^t 


f  (t)  =  sin(27Tfot)  exp  \  - 


where  f0  is  the  predominant  frequency  and  the  factor  o 


{ 


1 


i  3 1 
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controls  the  damping. 

In  order  to  determine  the  seismic  response  it  is 
necessary  to  evaluate  S(kj,z,t)  for  each  root  k,  over  the 
z-t  grid.  The  response  at  any  given  depth  z0  and  epicentral 
distance  r0  is  found  from  the  inverse  transform 


2.25 


00 


=  —  I 

a  1=1 


S  (k±,z,t)  J1  (k^) 
(J  (k.a))2 

O  1 


where  the  relation  J2(kja)  =  Ji(kja)  =  Jo(kta)  has  been  used 
(Abramowitz  and  Stegun  (1968)).  In  solving  for  S(kifz,t), 
equations  2.19  and  2.23  are  used  in  conjunction  with  the 
initial  conditions  imposed  by  2.13.  The  error  introduced  by 
truncation  of  the  series  is  of  the  order  2-3%  over  100X 
provided  a  sufficient  number  of  roots  k(  are  considered.  For 
the  case  f0=1  and  o= 4,  a  good  estimate  for  the  number  of 
terms  in  equation  2.25  is  "4a”.  For  the  P-SV  case,  extensive 
investigation  of  the  technique  places  the  number  of  terms  at 
"9a"  (Personal  communication  B.  G.  Mi khai lenko ) . 


2.5  Numerical  Considerations 

The  number  of  calculations  required  can  be  reduced  by  a 
consideration  of  the  rate  at  which  the  disturbance 
propagates  through  the  medium.  The  finite-difference 
equations  are  to  be  solved  for  all  depths  z=k(Az),  0^k<m,  at 
a  given  time  step.  The  points  in  z  to  which  the  disturbance 
has  not  yet  propagated  may  be  omitted  from  the  calculation 
in  the  initial  portion  of  the  grid  since  S(kifz,t)sO  in  that 


(  , 


. 
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region . 

In  addition,  for  a  given  receiver  depth  z0,  it  is 
unnecessary  to  compute  S(kj,z,t)  at  those  grid  points  that 
will  influence  the  disturbance  after  T,  the  length  of  the 
seismogram.  A  consideration  of  the  grid  elements  used  in  the 
finite-difference  calculation,  Figure  2.2,  shows  that  a  unit 
step  staircase  can  be  used  from  the  lower  grid  boundary  to 
the  depth  z0  in  excluding  grid  points  from  the  calculation 
in  the  latter  portion  of  the  grid  as  illustrated  in  Figure 
2.3. 


2.6  Stability 

Because  there  is  a  possibilty  of  unlimited 
amplification  of  errors  by  the  finite-difference  method  for 
arbitrary  choice  of  Az  and  At,  the  stability  criteria  for 
the  technique  must  be  determined.  The  von  Neumann  condition 
for  stability  can  be  applied  because  of  the  separability  of 
variables  in  the  problem.  In  this  method  a  harmonic 
decomposition  of  the  error  E  is  made  at  grid  points  at  a 
given  time  level.  The  form  of  the  error  function  E  can  be 
written  as  (Mitchell  (1977)): 

2.26 

i6  .  Az 

E  =  l  A.e  3 

j 

i6  .  (kAz) 

-  T  A.e  3 

j  3 


where,  in  general,  the  frequencies  0 j  are  arbitrary.  It  is 


. 


i-1  i  i+1 


z 


Figure  2.2  Grid  points  used  in  the  finite-difference 

calculation . 
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Figure  2.3  Grid  points  used  to  calculate  the  response. 
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necessary  to  consider  only  the  single  term  e  where 

0  is  any  real  number.  For  convenience,  suppose  that  the  time 
level  being  considered  corresponds  to  t=0.  To  investigate 
the  error  propagation  as  t  increases,  it  is  necessary  to 
find  a  solution  of  the  finite-difference  equation  which 
reduces  to  when  t  =  0.  Let  such  a  solution  be 


at  iBkAz 
e  e 

where  a=a(0 )  is,  in  general,  complex.  The  original  error 
will  not  grow  with  time  provided  that  iea^l<l  is  satisfied 
for  all  a.  To  determine  if  the  error  grows  with  time,  let 
the  error  be  written  in  the  form 
2.27 

_n  anAt  igkAz 
E,  =  e  e 

k 

.n  iBkAz 
=  A  e 


where  A=ea^t  .  The  quantity  A  is  referred  to  as  the 
amplification  factor  of  the  finite-difference  equation.  To 
facilitate  substitution  of  the  error  function  into  equation 
2.19  the  variable  coefficients  are  frozen  to  their  values  at 
a  specific  grid  point  and  the  stability  condition  determined 
from  the  von  Neumann  criterion  by  testing  the  modified 
equation  with  constant  coefficients.  The  test  should  in 
theory  be  carried  out  for  all  points  in  the  region  under 
consideration . 
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Writing  v2  =  m/p  and  a  =  M  equation  2.19  can  be  written 


as 

2e28 


,i+l 


•  -m  y  <i.i  -  H®2-2  *  *>’’> 


+ 


2/  At 
Vs(  Az 


Substitution  of  the  error  function  into  this  equation  yields 
2e29 


,n+l  iBkAz  _  2/ At\2  i6(k-l)Az 

vsVAz/ 


+  {2-2vs(§)2  -  k2v2(At)2}xneiekAz 


+  \n  i3  (k+1)  Az  _  ^n-lei6kAz 

vs\Az )  e 


2.30 


J 


{2  -  -  ki  vs<At>2} 

=  vs(H)2  {2cos6Az  -  2j  +  2  -  k2v2(At) 


2  +  2 


l/v2(At) 


2.31 


•H  M 
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Equation  2.31  may  be  written  in  the  form 


X  -  (X1+X2) X  +  XXX2  =  0 


where  and  X2  are  the  roots.  For  stability  it  is  necessary 
that  | X ! | < 1  and  | X 2 | < 1 ,  hence 


I  X !  +  X  2  |  *  |X1|  +|X2 


<  2 


Therefore 


r2{*±  - 


si  Az 


4sin2  Ml)  +  2 


2  2 

k^v  (At) 
is 


<  2 


«  „  2 /At \2  .2  2  ,  A  »  2 

2  -  4  v_  (  tt  )  -  k..  v_(At) 


s  \Az  ) 


l  s 


<  2 


The  corresponding  inequalities  are 


2.32 


vs(§)2  +Xvs(At>2  >  0 


2.33 


k2 

2/At\2  .  i  2 


Vs \ Az" )  +T~  vs(At)  <1 


The  stability  condition  is  given  by  the  last  inequality 
since  the  first  one  is  trivially  satisfied. 


. 


3 .  P-SV  CASE 


3.1  Equations  of  Motion 

The  equation  of  motion  for  a  vertically  inhomogeneous 
medium  is  given  by 
3.1 

••  ^ 
pu  =  (X+2y) V (V  *u)  +  ( V  *u) VX  -  y(VxVxu) 

+  (VyV)u  +  V(Vyu)  . 


Assuming  that  Lame's  parameters  are  functions  of  depth  only 
u-n{z)  ,  X  =  X(z),  and  (X  +  2/i )  =77  ( z )  the  terms  of  equation  3.1 
have  the  following  forms  in  cylindrical  coordinates  where 
u  =  (ur,u^,uz).  The  radial  and  vertical  components  of  the 
displacement  vector  u  are  given  by  ur  and  u^  respectively 
and  may  be  expressed  as 
Radial  Component: 


3.2 


n 


JW  1  3 

3r  V  r  cTr 


+ 


2  in 

3  z 


+ 


P 


3 


2 


u 


r 


Vertical  Component: 

3 . 3  3  /  ^z \  s  /n  i 

3z  y  n  3z  J  3z  \  r  3r 


(rur) 


n  JL 

r  3r 


r 


-  2 


in 

3  z 


1  3  (ru  ) X  u  3 
r  37  J  r  3r 


24 
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It  should  be  noted  that  in  the  original  wave  equation  a 
body  force  term  must  also  be  included.  For  the  case  of  a 
vertical  point  force  at  ( r  ,  z  )  =  ( 0  ,  d )  having  time  dependence 


f  (t) 
3.4 


? (r ,  z  ,  t) 


f  (t)  6 (r) 6 ( z-d)  ; 
2ur 


and  for  an  explosive  source 

3.5 


?<r,z,t)  =  f  (t)  7  (r2^rZ~~d>}- 


Equations  3.2  and  3.3  determine  the  displacement 
components  subject  to  initial  conditions 

3.6 


u 


3u 

at 


t=0 


=  0 


t=0 


and  boundary  conditions  at  the  free  surface 

3.7 


z=0 


3.8 


rz 


=  0 


z=0 


In  cylindrical  coordinates  the  boundary  conditions  at 
the  free  surface  can  be  written  as 


- 


.  , 

, 
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3.9 


a 


zz 


T7 


3.10 


As  for  the  SH  case,  the  dimensionality  of  the  problem 
can  be  reduced  through  the  introduction  of  finite  Hankel 
transforms.  The  choice  of  the  order  of  Hankel  transform  to 
be  applied  is  determined  largely  by  the  form  of  the  pair  of 
differential  equations  which  describe  the  motion,  and  the 
boundary  conditions  imposed  at  r=a. 

The  transforms  appropriate  to  the  boundary  value 
problem  expressed  by  3.2  and  3.3  with  boundary  conditions 
u  I  =0  and  3uz/3r|  =  0  involves  a  consideration  of  the 

Sturm-Liouvi lie  transforms  appropriate  to  this  differential 
system.  A  finite  first  order  Hankel  transform  is  applied  for 
ur  ,  and  a  finite  zero  order  transform  of  the  second  kind  is 
applied  to  uz. 

The  transforms  are  given  by 

3.11  ra 

S(ki#z,t)  =  ur  (r,z ,t) J1 (kir) rdr 

o 


3.12 


R (ki , z ,t) 


-a 

u  (r,z,t) J  (k.r) rdr 
z  o  1 

o 


. 
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where  kj  are  the  roots  of  J1(kia)=0  (Sneddon  (1972)). 

Direct  application  of  this  pair  of  transforms  to 
equations  3.2  and  3.3  reduces  all  operations  with  the 
variable  "r"  to  scalar  multiplication  in  the  transform 
domain.  Thus  in  P-SV  problems  with  a  point  source  on  r  =  0,  a 
zero  order  Hankel  transform  is  appropriate  for  uz  and  ozz 
whereas  a  first  order  transform  is  required  for  uf  and  t  . 
(Aki  and  Richards  Vol  I  (1980)). 

The  inverse  transform  for  the  radial  component  is 

3.13 

®  S  (k.  ,  z  ,  t)  J,  (k .  r) 

u  (r ,  z , t)  =  4  l  - ± - 4^— 

r  a2  i=l  (J2(kia)) 


and  for  the  vertical  component 


3.  14 


«  R(k  .  ,  z  ,t)  J  (k  .  r) 

u  (r ,  z,t)  =  4  l  — i - 

2  a2  i=l  (J  (k.a)) 


The  inversion  for  the  radial  component  can  be  expressed  in  a 
form  similar  to  that  for  the  vertical  component 


3.15 


u  (r,z,t)  =  -y  l 
r  a2  i=2 


2  »  S  (k^fz,t) (k^r) 


(J  (k.a) ) 
o  1 


where  the  relation 

-J0 (k.a)  =  J  (k.a) 

2  i  o  l 


(J1(kia) 


has  been  used  (Abramowitz  and  Stegun  (1968)). 
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As  for  the  SH  case  the  introduction  of  boundary 
conditions  at  r=a  introduces  a  fictitious  reflecting 
surface.  The  boundary  r=a  must  be  chosen  sufficiently  far 
that  reflections  from  it  do  not  influence  the  seismic 
response  in  the  region  of  interest. 

The  transformation  of  each  term  in  the  equation  for  the 
radial  component  involves  integration  by  parts: 


a 


o 


>a 

4 

o 


•a 

4 

o 


^1  (^ir )  rdr  =  R(ki,z,t) 

J-(k.r)rdr  =  -k.  -3-  (yR(k.,z,t)) 
±1  l  o  z  i 

(k±r ) rdr  =  -ki  |£  R(ki,z,t) 


•a 

4 

o 


1/11 
n  Sr  \r  Sr 


J1(kir)rdr 


-nk?  S(ki 


Z,t)  . 


The  equation  for  the  radial  component  becomes 
3.16 

4  0  If)- H  +  ki  h  (yR>  -  2ki  H  R-nkis  - 


s2s 

st2 


To  transform  the  equation  for  the  vertical  component, 
consider  the  terms  of  equation  3.3. 


3  ^(r  Jr  (rur})  Jo(kir)rdr  -  k±  Bz  (nS(ki,z,t) 


o 
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o 


V  3 
r  Tr 


“a  3 Vi  1  3 


Jz  r  "cTr 


o 


o 


ra 


p  3 
r  Jr 


The  equation  for  the  vertical  component  becomes  upon 
transformation : 

3.17 


2  * 


3 1 


To  determine  the  radial  and  vertical  displacement 
components  it  is  necessary  to  solve  the  coupled  system  given 
by  equations  3.16  and  3.17  using  finite-difference  methods. 
The  inverse  transforms  given  by  3.14  and  3.15  define  the 
components  ur  and  u^  respectively. 

3.2  Finite-Difference  Equations 

It  remains  to  develop  finite-difference  analogues  for 
equations  3.16  and  3.17.  In  the  equation  for  the  radial 
component  the  terms  in  R(kj,z,t)  can  be  written  in  the  form 
( equat i on  3.16): 


3.18 
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ik  (<rWr,k-i-pk+i-uk-i)  (Rk+rRk-i) 

+  ^k+r^k-l’  ^k+l^k-l’} 


k. 

_ i_ 

2Az 


(,1k-2l,k)+uk+i 


R 


k+1 


(r>k-2yk)+yk-D 


k-1, 


Arithmetic  averages  have  been  used  in  order  to  provide 
some  smoothing  of  the  medium  parameters  and  response.  In  a 
similar  manner,  the  terms  in  equation  3.17  for  the  vertical 
component  which  involve  S(kifz,t)  can  be  written: 

3.19 

u  1  3  _  ..  3S  _  ,  _ 

dZ 


<ns)  -  y  |f  -  2  U  s 


'  (  — 
i  \  9z 

=  kx  { H  +  k  (n_2p)s} 


=  h.  f 

4  Az  ^ 


(rlk+l+rlk-l_,Jk+l'wk-l>  <Sk+l  Sk-1) 


+  (nk+l‘2uk+l"rik-l+2l,k+l)  ^k+l^k-l1 


Al 

2  Az 


(rik+r2pk+i)+l,k 


k+1 


(nk-r2yk-i)+yk 


k-1 


In  the  same  manner  as  for  the  SH  case,  the  following 
analogues  are  used: 

3.20  9  /  3S\  akSk-l  ”  (ak+ak+l)Sk  +  ak+lSk+l 

^(2)  97  >  = 


87 


(Az) 


I 
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where  a, 
k 


2tv 

Pk+y 


k-l 

k-1 


3.21 


Jz  '  n (z) 


5R 

dZ 


hk^-l  ”  (bk+bk+l 


)Rk  +  bk+lRk+l 


(Az) 


where  b. 


2nknk-l 

nk+nk-l 


3.22 


-2sk+sk'1\ 
(At)2  / 


.i+1 


3.23 


dS 

Y  Ft  = 


-  S 


i-1 


Y, 


2  (At) 


The  analogues  for  the  time  derivative  in  R  are 
indentical  to  3.22  and  3.23,.  Making  use  of  these  relations 
the  finite-difference  equivalent  of  equation  3.16  for  the 
radial  component  is: 

3.24  a,  si  ,  -  (a, +a,  .  )  S.1  +  a.  ..si,,  /  si+1-2si+si_1 

k  k-l  k  k+1  k  k+1  k+1 


(Az) 


~  P, 


’k  *k 


(At) 


-k2n,  si  - 


‘i  k  k  2  Az 


2 Az  ((\~2yk+Pk+l)Rk+l  “  (nk"2pk+yk-l)Rk-l}) 


-  y  =0 

Tk  \  2 At  1  * 


* 
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This  equation  is  to  be  solved  for  S**1  explicitly  and 
can  be  written  as 
3.25 


Yk  \  „i+l  _  akSk+l 


(At) 


2  2  At  k 


(Az) 


!*2*±i  +  k^  .  JL 

2  1  (Az)2  1  k  (At) 


_ )  s'1 

2  bk 


a.sr  ,  k.  • 

A  *9  (hi,“2y, +p,  .  )Ri. 


(Az) 


2  2Az  k  k  k+1  k+1 


k. 


+  2  Az  (T1k"2yk+yk-l)Rk-l 


(At) 


_  _  2lL)  s1’1 

2  2  At  k 


l<k^m 


The  finite-difference  analogue  of  the  equation  for  the 
vertical  component  is 
3.26 

bkRk-i  -  (bk+bk+i>4 +  bk+iRk+i  .  ( Rki+1-2R^+Rki_1 
- 2 - pk  - JTT ^ - 


(Az) 


(At) 


“kiykRk  +  2  Az  (nk+l~2yk+l+yk)Sk+l“(rik-l~2yk-l+yk)Sk-l 
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or  equivalently  as 
3.27 
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pk 

(At) 


2 


+ 


b 


R1 

k+1  k+1 
(Az)2 


/VVl 

\  (Az)2 
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,2 
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i  k 


+  ,2  +  2  Az  (T1k+l“2yk+l+pk)  Sk+1 
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Equations  3.25  and  3.27  are  used  to  evaluate  successive 
values  of  R(kifz,t)  and  S(kj,z,t).  The  displacements  ur  and 
uz  are  found  from  the  inverse  transforms  given  by  equations 
3.14  and  3.15.  It  should  be  noted  that  the  body  force  terms 
must  also  be  incorporated  into  the  finite-difference 
equations . 


3.3  Source  Terms 

For  the  case  of  a  vertical  point  force,  the  transform 
of  the  body  force  term  is  included  on  the  left  side  of 
equation  3.17  for  the  vertical  component.  The  transform  is 
given  by 
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3.28 


F  (r,z ,t) JQ (k^r ) rdr 
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f  (t)5  (z-d) 

2tt 


6  (r) JQ (kir) dr 


o 


f  (t)6(z-d) 

4tt 


This  term  is  to  be  added  to  R(kifz,t)  at  depth  z=d.  The 
Dirac  delta  function  in  this  case  has  been  approximated  by 
1/Az  since  6(z-d)*Az=1.  The  function  f(t)  is  an 
exponentially  damped  sine  function  as  for  the  SH  case. 

For  an  explosive  source 


3.29 


F  (r ,  z  ,  t)  =  f  (t)  V 


=  f(t) 


_3_ 

dz 


6  (r) 6 (z-d) 
2iTr 


The  transform  of  the  radial  component  is 


3.30 


f(t> 


6 (r) 5 (z-d) 
2iTr 


)  J1 (kir) rdr 
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kif  (t)5(z-d) 


4  TT  Az 


This  term  is  to  be  added  to  S(kj,z,t)  at  depth  z=d.  The 


N  > 


.  1 


transform  of  the  vertical  component  for  this  case  is 
3.31 
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f  (t) 


_3_  /  6  (r)  6  (z-d) 
8z  \  2mr 


J  (k.r)rdr 
o  1 


=  £  «  (Z-d) 


f  (t) 
4  TT 


If  5(z-d)  is  approximated  by  1/Az  at  depth  z  =  kAz  then 
an  approximation  to  d/dz  6(z-d)  is  given  by  1/(Az)2  at 
z  = ( k  —  1 ) Az  and  -1/(Az)2  at  z=(k+1)Az.  The  terms  to  be  added 
to  R ( k  j , z , t )  are 
3.32 

— f  at  z  =  (k-l)Az 

4 7T  ( Az ) 


f  (t) 
4tt(Az)2 


at  z  =  (k+1) Az 


for  the  duration  of  the  seismic  pulse. 


3.4  Free  Surface  Condition 

The  finite-difference  equations  given  by  3.25  and  3.27, 
with  the  incorporation  of  the  appropriate  body  force  terms, 
are  to  be  solved  for  all  grid  steps  k  where  1<k<m.  At  the 
free  surface,  k=0,  the  assumption  of  a  homogeneous  medium 
allows  the  equations  of  motion  to  be  written  in  a  simplified 
form.  Homogeneity  for  k=0,1,2  enables  equations  3.16  and 
3.17  to  be  written  as 
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3.33 


3.34 


where  v2  =  p/p  and  v2  =  (X+2 u)  /  p ,  If  finite-difference 

S  p 

analogues  are  written  for  these  equations,  they  will  contain 
the  fictitious  values  Ri ,  and  S I i .  In  order  to  eliminate 
these  values  from  the  equations,  the  boundary  conditions  at 
the  free  surface  given  by  equations  3.7  and  3.10  may  be 
used.  From  equations  3.7  and  3.9  we  may  write,  for  ozz 


Application  of  a  zero  order  finite  Hankel  transform  yields 
for  z=0: 
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a  3u 
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n  J  (k  .  r)  rdr  =  0 

3z  o  l 
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XkiS(ki,z,t)  +  n  gy  (ki,z,t)  =  0 
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The  finite-difference  analogue  to  this  last  equation 
yields 
3.36 


In  a  similar  manner,  the  equation  for  the  shear  stress  at 
the  free  surface  given  by  3.8  and  3.10  can  be  used  to 
determine  a  substitution  for  the  fictitious  value  S 1 i . 
Application  of  the  transform  yields 
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The  finite-difference  analogue  to  this  free  surface  stress 

condition  is 

3.38 

S1,  =  -  2k • AzR^ . 

-11  l  o 


From  3.33  the  finite-difference  relation  for  z=0  may  be 
obtained  as 
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The  explicit  solution  for  Si*1  becomes 
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Substitution  of  equation  3.38  to  eliminate  S i w  we  have 
3.39 
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The  explicit  solution  for  Ri*1  is  obtained  from 
equation  3.34  in  the  same  way: 
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The  explicit  solution  for  Ri*1 
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Substitution  of  equation  3.36  for  Ri ,  yields 
3.40 
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Making  use  of  the  relations  given  by  3.39  and  3.40,  the 
truncation  error  at  the  free  surface  as  well  as  throughout 
the  interior  of  the  grid  is  of  the  order  0 ( ( Az ) 2 , ( At ) 2 ) .  The 
truncation  error  at  the  free  surface  was  thought  to  be 
responsible  for  the  generation  of  the  S*  non-geometric 
arrival  identified  by  Hron  and  Mikhailenko  (1981).  A  Taylor 
series  expansion  at  the  free  surface  for  R(kifz,t)  and 
S(kj,z,t)  was  then  used  to  reduce  the  truncation  error  to 
the  order  of  approximation  0((Az)\  (At)2).  The 

incorporation  of  this  improved  boundary  condition  did  not 
alter  the  results  substantially  and  the  nature  of  the 
arrivals  was  unchanged.  The  basic  features  of  the  S*  arrival 
were  first  inferred  by  Hron  and  Mikhailenko  (1981)  from 
synthetic  seismograms  based  on  the  Alekseev-Mi kha i lenko 


method . 
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3.5  Stability 

In  order  to  determine  the  stability  condition  for  the 
P-SV  case  the  coupled  finite-difference  equations  to  be 
analyzed  by  the  von  Neumann  method  are  given  by  equations 
3.25  and  3.27. 

In  the  same  manner  as  for  the  SH  case,  if  we  assume 
that  Yk=0  and  set  the  variable  coefficients  to  their  values 
at  a  specific  grid  point,  the  modified  equations  in  terms  of 
shear  and  compressional  velocities  are 
3.41 
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where  v2  =  \i/ p  and  v2  =  77/p . 

a  p 

For  a  system  of  coupled  equations,  Richtmyer  and  Morton 
(1967)  dicuss  in  detail  the  basis  for  using  a  harmonic 
decomposition  in  terms  of  a  Fourier  series.  In  our 
particular  case  it  will  be  necessary  to  determine  the 
amplification  matrix  from  which  the  characteristic  equation 


. 


yields  the  eigenvalues.  To  ensure  stability  the  spectral 
radius,  or  maximum  of  | X  s  |  i=1,2,3,4  must  be  uniformly 
bounded :  |  X  ;  |  <  1  . 

Given  a  properly  posed  initial  value  problem  and  a 
finite-difference  approximation  to  it,  Lax's  Equivalence 
Theorem  states  that  stability  is  the  necessary  and 
sufficient  condition  for  convergence.  The  consistency 
condition  requires  that  the  truncation  error  **0  as  (At)->0 
for  O^t^T.  The  amplification  factors  or  eigenvalues  can  be 
obtained  by  substituting  the  error  function  directly  into 
the  difference  equations.  This  yields  a  system  of  linear 
equations  for  Rw  *  1  and  S**1  and  the  vanishing  of  their 
determinant  determines  X.  In  our  case,  the  coupled  equations 
given  by  3.41  and  3.42  can  be  expressed  as 
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where  Pk+1  =  Si  and  Qk*1  =  Rj  have  been  used  to  reduce  the 
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system  to  a  two-level  scheme.  If  we  write 
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then  the  set  of  equations  can  be  written  in  the  form: 
3.45 
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In  terms  of  an  amplification  matrix,  W=(P,Q,R,S)  can  be 


expressed  as: 
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3.46 


where 


/  v  At\2  9  2  2  2 

a  =  |  ■  (-4 sin  ( 3 Az/2 ) )  +  2  -  k^vs ( At) 

ik  .  (At) 2  o  ? 
b  =  — ^ - (v  V)  sin  (  3Az) 


3.47 


c  =  -b 


v  At\2  2 

d  =1  I  ( -4sin  ( 3 Az/2 ) ) 


+  2  - 


.  2  2  /  a  i  \  2 

kivp(At) 


The  eigenvalues  of  the  amplification  matrix  in  equation 
3.46  are  given  by  | A— X I | =  0 .  The  amplification  factors  are 
found  from  the  roots  of 

3.48  .  9  2 

A  -  (a+d)A  +  (2+ad-bc) A  -  (a+d) A +  1  =  0. 

If  we  denote  the  eigenvalues  as  { X  ,  ,  A  2  ,  A  3  ,  A  „ }  then  the 
equation  of  degree  4  having  these  roots  can  be  written 
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Comparison  of  3.48  and  3.49  yields  the  relation 
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The  requirement  for  stability  is  | X  j  | < 1 .  Imposing  this 

condition  on  the  eigenvalues  we  have 
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From  equation  3.50  the  requirement  for  stability  is 
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The  two  equivalent  inequalities  are 
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The  first  inequality  is  trivially  satisfied  since  all  terms 
are  positive  definite.  The  necessary  and  sufficient 
condition  for  stabilty  of  the  finite-difference  scheme  is 
given  by  the  last  inequality. 


4.  NUMERICAL  MODELS 


In  order  to  provide  a  comparison  of  results  with  those 
presented  in  the  literature,  a  particular  model  has  been 
duplicated  for  the  case  of  SH  waves  and  then  extended  for 
other  source  types.  The  model  used  for  the  calculation  of 
synthetic  seismograms  for  a  vertically  inhomogeneous  medium 
is  one  presented  by  M.  Korn  and  G.  Muller  (1983)  in  which 
two  thin  low-velocity  layers  representing  coal  seams  are 
embedded  in  a  homogeneous  half-space.  In  their  paper  SH 
waves  are  generated  by  a  horizontal  single  force  at  the 
surface  and  the  results  confirmed  by  comparison  with  the  ray 
reflectivity  method. 

The  statement  of  the  problem  is  given  in  cylindrical 
coordinates  (r,0,z)  with  the  z-axis  pointing  downward  into 
the  vertically  inhomogeneous  half-space  whose  surface  is  at 
z=0.  As  before  the  material  parameters  are  the  density  p( z), 
rigidity  u( z)  and  S  velocity  P ( z )  .  At  the  surface  a 
horizontal  single  force  f(t)  acts  in  the  direction  0=0.  The 
boundary  values  of  the  tangential  stresses  are 
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The  differential  equation  describing  u^,  the  azimuthal 
displacement  is 
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Figure  4.1  Velocity  vs  depth  for  the  coal  seam  model. 
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The  stress  relation  tzq  =  p(3u^/3z)  together  with  4.2  yields 
4.4 
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Since  the  azimuth  <t>  does  not  appear  in  4.3,  equation  4.4 
implies  u^  ~  sin  <p  for  constant  r  and  z.  Hence 
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As  before,  the  initial-boundary-value  problem  to  be 
solved  is  given  by 
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The  finite  Hankel  transform  pair  used  in  reducing  the 
dimensionality  of  the  problem  is  given  by 


4.9 


S  (ki  ,  z  ,  t ) 


■a 

u(r,z,t) JQ (k^r) rdr 


o 


49 


4.10 
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where  kj  are  the  roots  of  Jo(kia)=0.  The  differential 
equation  for  the  Hankel  transform  S(kj,z,t)  and  the  boundary 
and  initial  conditions  are  given  by 
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Upon  solution  of  this  transformed 
initial-boundary-value  problem,  the  displacement  is  given  by 
the  inverse  transform  4.10. 

For  the  numerical  calculations,  the  coal  seam  model 
involves  layers  each  2m  thick  having  density  1.6  g/cm3  and 
S-velocity  866  m/s  embedded  in  a  half-space  of  density  2.6 
g/cm3  and  S-velocity  1732  m/s.  The  coal  seams  are  placed  at 


50 


depths  of  200m  and  250m.  On  a  profile  at  the  surface,  along 
<p=90°,  one  expects  direct  waves,  primary  seam  reflections, 
interseam  multiples,  and  multiples  between  the  surface  and 
the  seams. 

In  order  to  match  the  arrivals  given  in  their  paper, 
the  dominant  frequency  of  the  force  f(t)  was  taken  to  be  60 
Hz  rather  than  30  Hz  as  stated  by  Korn  and  Muller.  In  the 
seismograms  presented  here,  epicentral  distances  and  depths 
are  given  in  terms  of  S-wavelengths  for  the  half-space.  In 
all  figures,  times  are  given  in  periods  rather  than  seconds. 
With  these  changes  the  coal  seams  are  at  depths  of  6.90  WL 
and  8.625  WL.  One  set  of  10  receivers  is  located  on  the 
surface  of  the  half-space  at  epicentral  increments  of  0.862 
WL.  For  the  vertical  seismic  profiles  another  set  of  24 
receivers  at  equal  depth  increments  of  0.375  WL  has  been 
placed  at  an  epicentral  distance  of  4.31  WL.  For  the  P-SV 
case  the  source  is  located  at  a  depth  of  0.75  WL.  Sources 
for  the  SH  case  are  invariably  located  on  the  surface. 

4.1  Horizontal  Point  Force  -  Surface  Traces 

In  Figure  4.3,  the  direct  arrival  A  as  well  as  the 
primary  reflections  B  and  C  are  evident.  In  order  to 
identify  each  of  the  remaining  arrivals,  the  direct  arrival 
has  been  windowed  out  and  the  later  arrivals  amplified  as 
shown  in  Figure  4.4.  The  primary  reflections  B  and  C, 
interseam  multiples  D  and  E,  and  the  secondary  reflections  F 
and  G  have  been  identified  and  confirmed  by  matching  travel 
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Figure  4.2  Traces  from  Korn  and  Muller. 
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Figure  4.3  Horizontal  Point  Force  -  Surface  Traces 
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Figure  4.4  Horizontal  Point  Force 
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times  using  elementary  ray  methods. 

There  is  a  noticeable  distortion  of  the  original  pulse 
shape  in  the  primary  reflections  and  later  arrivals.  This  is 
due  to  the  complex  interference  phenomena  present  in  the 
wavelets  which  are  reflected  from  or  transmitted  through  the 
very  thin  layers  which  represent  the  coal  seams. 

4.2  SH  Point  Torque  Source  -  Surface  Traces 

Only  the  windowed  seismograms  are  presented  in  Figure 
4.5  for  this  case.  The  significant  difference  between  these 
arrivals  and  those  for  the  horizontal  point  force  is  the 
amplitude  dependence  of  the  arrival  on  the  angle  e  where  © 
is  the  angle  made  by  the  ray  with  the  z-axis.  Ben-Menahem 
(1981)  confirms  that  the  amplitude  of  the  disturbance  varies 
directly  with  sin  ©  for  a  point  torque  source. 

4.3  SH  Point  Sources  -  Vertical  Profiles 

From  the  vertical  profiles  given  in  Figures  4.6  and 
4.7,  the  nature  of  the  wave  propagation  is  evident.  Energy 
may  arrive  at  a  particular  depth  either  from  above  or  below. 
The  direct  arrival  A  represents  the  downgoing  wavefield; 
arrivals  B  and  C  are  reflections  from  each  of  the  seams  and 
represent  upward  travelling  disturbances.  The  interseam 
multiples  are  a  result  of  the  reflection  labeled  D.  The 
primary  reflection  B  from  the  uppermost  coal  seam  is 
reflected  once  again  at  the  free  surface  and  propagates 
downward  into  the  half-space  as  arrival  E. 
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Figure  4 
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Figure  4.6  Horizontal  Point  Force 
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Figure  4.7  SH  Torque  Point  Source 
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4.4  Explosive  Source  -  Surface  Traces 

The  set  of  seismograms  at  the  free  surface  are  given  in 
Figures  4.8  and  4.9.  The  direct  arrival  A  as  well  as  the 
reflections  from  the  seams  given  by  C  and  D  are  prominent  in 
Figure  4.8.  The  weak  arrival  B  is  a  low  amplitude  Rayleigh 
wave  as  confirmed  by  its  apparent  velocity. 

4.5  Explosive  Source  -  Vertical  Profiles 

The  vertical  profiles  as  given  in  Figures  4.10  and  4.11 
illustrate  clearly  the  propagation  of  the  disturbance  in  the 
half-space.  The  direct  P  arrival  identified  as  A  undergoes 
significant  conversion  to  SV  upon  reflection  at  the  free 
surface.  The  resulting  S  wave  is  identified  by  B.  Arrivals  C 
and  D  are  reflections  from  the  first  coal  seam,  their 
arrival  times  -at  the  free  surface  matching  those  given  in 
Figure  4.8  which  were  again  confirmed  by  ray  methods. 

4.6  Vertical  Point  Force  -  Surface  Traces 

The  direct  P  arrival  as  well  as  a  very  strong  Rayleigh 
wave  are  the  most  discernable  features  in  Figures  4.12  and 
4.13.  For  the  vertical  component  in  particular,  as  given  in 
Figure  4.13,  the  Rayleigh  wave  is  the  outstanding  feature. 
The  nature  of  the  radiation  pattern  from  the  vertical  point 
force  is  responsible  in  this  case  for  its  strength.  The 
force  acts  as  a  source  of  both  P  and  S  waves  and  the  strong 
vertical  component  enhances  the  magnitude  of  the  surface 
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Figure  4.8  Explosive  Source  -  Horizontal  Component  of 

Surface  Traces 
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Figure  4 


.9  Explosive  Source  -  Vertical  Component  of  Surface 

Traces 
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Figure  4.10  Explosive  Source  -  Horizontal  Component  of 

Vertical  Profile 
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Figure  4.11  Explosive  Source  -  Vertical  Component  of 

Vertical  Profile 
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Figure  4.12  Vertical  Point  Force  -  Horizontal  Component  of 

Surface  Traces 
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Figure  4.13  Vertical  Point  Force  -  Vertical  Component  of 

Surface  Traces 
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4.7  Vertical  Point  Force  -  Vertical  Profiles 

Two  features  are  immediately  evident  in  Figures  4.14 
and  4.15.  The  first  is  that  there  is  both  P  and  S  radiation 
from  the  source.  This  is  confirmed  by  the  separation  of  the 
direct  P  and  S  arrivals  on  the  surface  trace.  The  second 
feature  is  the  more  complex  nature  of  the  wavefield  produced 
by  this  source  than  by  the  explosive  source  as  a  result  of 
both  P  and  S  waves  emanating  from  the  source  and  partial 
conversions  at  the  free  surface  upon  reflection. 

4.8  Conclusion 

For  the  generation  of  synthetic  seismograms  the 
Alekseev-Mi khailenko  method  has  the  advantage  of  providing  a 
seemingly  true  representation  of  the  wavefield  produced  by 
various  point  sources.  Included  in  the  traces  are  found 
direct  arrivals,  reflections  and  their  multiples,  and 
surface  waves.  An  additional  advantage  is  the  numerical 
accuracy  provided  through  the  use  of  Hankel  transforms  in 
reducing  the  dimensionality  of  the  wave  equation  thereby 
limiting  grid  dispersion.  The  use  of  this  numerical  method 
in  vertical  seismic  profile  studies  is  evident,  however  ray 
methods  must  still  be  used  to  identify  particular  arrivals. 
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Figure  4.14  Vertical  Point  Force  -  Horizontal  Component  of 

Vertical  Profile 
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Figure  4.15  Vertical  Point  Force  -  Vertical  Component  of 

Vertical  Pro'file 
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